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Cn ' Abstract. A new computation method of frequentist p-values and Bayesian 

posterior probabilities based on the bootstrap probability is discussed for the 
multivariate normal model with unknown expectation parameter vector. The 
null hypothesis is represented as an arbitrary-shaped region. We introduce 
new parametric models for the scaling-law of bootstrap probability so that the 
multiscale bootstrap method, which was designed for one-sided test, can also 
computes confidence measures of two-sided test, extending applicability to a 

Mm ' wider class of hypotheses. Parameter estimation is improved by the two-step 

multiscale bootstrap and also by including higher-order terms. Model selection 
is important not only as a motivating application of our method, but also as 
an essential ingredient in the method. A compromise between frequentist and 
C^ ' Bayesian is attempted by showing that the Bayesian posterior probability with 

j-»s ' an noninformative prior is interpreted as a frequentist p-value of "zero-sided" 

test. 



1. Introduction 



> 

m 
in 

^^ I Let Y = (Yi, . . . , Y„i^i) be a random vector of dimension m + 1 for some in- 

teger m > 1, and y — {yi, . . . ,ym+i) G K™"*"^ be its observed value. Our ar- 

C^^ I gument is based on the multivariate normal model with unknown mean vector 

/x = (/xi, . . . , Hm+i) & M"'"'"^ and covariance identity Im+i, 

O' (1) i"-iV„,+i(/i,/^+i), 

where the probability with respect to ([T]) will be denoted as P(-|/i). Let Hq C M™+^ 

be an arbitrary-shaped region. The subject of this paper is to compute measures 

of confidence for testing the null hypothesis fi G Hg. Observing y, we compute 

C^ I a frequentist p- value, denoted p{HQ\y), and also a Bayesian posterior probability 

7r(ffol?/) with a noninformative prior density 7r(/i) of ii. 

Th is is th e pr oblem of region s discussed in literature: lEfron et all ( 19961 ). lEfron and Tibshirani 
(|l998l ). and [Shimodairai (|2002i [2004. 20081 . The confidence measures were calcu- 



lated by the bootstrap methods for complicated application problems such as the 
variable selection of regression analysis and phylogenetic tree selection of molecu- 
lar evolution. These model selection problems are motivating applications for the 
issues discussed in this paper, and the normal model of ([T]) is a simplification of 
reality. Let X = {xi, . . . , Xn} be a sample of size n in application problems. We 
assume there exists a transformation, depending on n, from A" to y so that Y is 
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approximately normalized. We assume only the existence of such a transformation, 
and do not have to consider its details. Since we work only on the transformed vari- 
able Y in this paper for developing the theory, readers may refer to the literature 
above for the examples of applications. Before the problem formulation is given in 
Section [21 our methodology is illustrated in simple examples below in this section. 
The simplest example of Hq would be the half space of IR™+^, 

(2) H^ : Mm+i < 0, 

where the notation Hq, instead of Hq, is used to distinguish this case from an- 
other example given in ([3]). Only Hm+i is involved in this Hq, and one-dimensional 
normal model Y„i+i ^ -/V(/i„j+i, 1) is considered. Taking iim+i > as an alter- 
native hypothesis and denoting the cumulative distribution function of the stan- 
dard normal as <&(•) with density </>(•), the unbiased frequentist p- value is given as 
PiHoly) = $(-y„,+i). 

A slightly complex example of Hq is 

(3) Hq : -d< ^im+i < 

for d > 0. The rejection regions are Hm+i > c and ym+i < ~d — c with a critical 
constant c, which is obtained as a solution of the equation 

(4) $(-c)-f $(-d-c) = a 

for a specified significance level < a < 1. The left hand side of Q is the rejection 
probability P{Ym^i > cV Ym+i < —d — c\fi) when fi is on the boundary of Hq, i.e., 
fjLm+i = or /im-i-i = —d. The frequentist p- value is defined as the infimum of a such 
that Hq can be rejected. This becomes p[HQ\y) = ^[—ym+i) + ^[—d — ym+i) for 
yrn+i > -d/2siiidp{Ho\y) = $(j/„,+i) + $(d+y,„+i) for y„j+i < -d/2. Considering 
the case, say, 

(5) d=l, 2/,„+i = -0.1, 

we obtain piH[,\y) = 0.540 and p{HQ\y) = 0.724. 

These two simple cases of Hq and Hq exhibit what lEfron and Tibshiranil (jl998| ) 
called paradox of frequentist p-values. Our simple examples of ^ and ^ suffice for 
this purpose, alth ough they had actually used the spherical shell example explained 



later in Section ID lEfron and Tibshiranil (J199S) indicated that a confidence measure 



should be monotonically increasing in the order of set inclusion of the hypothesis. 
Noting Hq C Hq, therefore, it should be p{HQ\y) < p{HQ\y), but it is not. This kind 
of "paradox" cannot occur with Bayesian methods, and 7r(i/o|y) < i^iH^^y) holds 
always. Considering the flat prior 7r(/i) = const, say, the posterior distribution of 
/i given y becomes 

(6) p.\y r^ N„^+i{y,I^+i), 

and the posterior probabilities for the case dS]) are 7r(iJo|y) = ^{—ym+i) = 0.540 
and 7r(i?o|2/) = ^(— 2/m+i) — $(— d— Vm+i) = 0.356. The "paradox" of frequentist 
p- values may be nothing surprise for a frequentist statistician, but a natural conse- 
quence of the fact that p(7Jq |y) is for a one-sided test and p{HQ\y) is for a two-sided 
test; The power of testing is higher, i.e., p- values are smaller, for an appropriately 
formulated one-sided test than a two-sided test. In this paper, we do not intend 
to argue the philosophical question of whether to be frequentist or to be Bayesian, 
but discuss only computation of these two confidence measures. 
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Computatio n of the confidence measures is made by the bootstrap resamphng 
of lEfronI ( 19791 ). Let X* — {x^, . . . , x*/} be a bootstrap sample of size n' obtained 



by resampling with re placement f r om X . The idea of bootstrap probability, which 
is introduced first by iFelsensteinI ( 1985[ ) to phylogenetic inference, is to generate 



X* many times, say B, and count the frequency C that a hypothesis of interest 
is supported by the bootstrap samples. The bootstrap probability is computed as 
C/B. Recalling the transformation to get y from X, we get Y* by applying the same 
transformation to X* . For typical problems, the variance of Y* is approximately 
proportional to the factor 



9 IL 



as mentioned in IShimodairal ( 20081 ). Although we generate X* in practice, we only 



work on y* in this paper. More specifically, we formally consider the parametric 
bootstrap 

(7) Y*\y^N^+iiy,a^Irn+i), 

which is analogous to ^ but the scale cr is introduced for multiscale bootstrap. 
The bootstrap probability is defined as 

(8) a^2{Ho\y) ^ P,2iY* e Ho\y), 

where P„2{-\y) denotes the probability with respect to ([7]). For computing a crude 
confidence measure, we set cr = 1, or n' = n in terms of X* , so that the distribution 
d?]) for Y* is equivalent to the posterior ^ for fi. This gives an interpretation 
of the bootstrap probability that ai (Ho\y) = Trfffnlw) for any ffp under the flat 
prior. In the multiscale bootstrap of IShimodairal ( 2002I . I2004L 120081 ) . however, we 



may intentionally alter the scale from cr = 1, or to change n' from n in terms of 
X* for computing p{Ho\y). Let ai, . . . , (Tm be M different values of scale, which we 
specify in advance. In our numerical examples, Af = 13 scales are equally spaced 
in log-scale between ai = 1/3 and ai^ = 3. For each i = 1, . . . , M, we generate X* 
with scale ai for Bi times, and observe the frequency Q. The observed bootstrap 
probability is a^2 = d/ Bi. 

How can we use the observed a^2, . . . ,a„2 for computing p(iJo|2/)? Let us 
assume that H^ can be expressed as ([3]) but we are unable to observe the values 
of y„i+i and d. Nevertheless, by fitting the model aa2{Ho\y) — $(— ym+i/f) — 
$(— (d + j/m+i)/'?') to the observed d^-s, . . . , a.^2 , we may compute an estimate (^ of 
the parameter vector ip — (j/m+i, d) with constraints d > and j/m+i > —d/2. The 
confidence measures are then computed as p(i?o|y) = *&(— j/m+i) + *&(— d — Vm+i) 
and 7r(_ffo|y) = ^{—Vm+i) — ^{—d— ym+i) ■ In case we are not sure which of ([2]) and 
dHI is the reality, we may also fit a„2{H'^\y) = $(— ym+i/cr) to the observed d^2's 
and compare the AIC values ( Akaikd . Il974h for model selection. In practice, we 



prepare collection of such models describing the scaling-law of bootstrap probability, 
and choose the model which minimizes the AIC value. 

2. Formulation of the problem 

The examples in Section [T] were very simple because the boundary surfaces of 
the regions are flat. In the following sections, we work on generalizations of ([2]) and 
([3]) by allowing curved boundary surfaces. For convenience, we denote y = (u, v) 
with u — {yi, . . . ,ym) and v — ym+i- Similarly, we denote /x = (0,/Zm+i) with 
9 = {fii, . . . ,iJ.m) G M™. As shown in Fig. [H we consider the region of the form 
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Ho = {{0,^^rn+l) I -d-h2{e) < n,n+i < -hi{e), 9 e M™}, where hi{e) and /12(e) 
are arbitrary functions of 6. This region will reduce to ^ if hi{9) = h2{6) = for 
all 9. The region may be abbreviated as 

(9) Ho:~d~ h2{9) < Ai™+i < -hi{9). 

Two other regions Hi : fim+i > ^hi{9) and H2 : fJ-m+i < —d— /12(e) as well as two 
boundary surfaces dHi : ^im+i = —hi{9) and dH2 : Hm+i = —d~ /12(e) are also 
shown in Fig. [1] We define Hq = Hq U H2, or equivalently as 

(10) i/^:Mm+i<-/ii(e). 

The boundary surfaces of the hypotheses are OHq — dHi U dH2 for the region Hq, 
and OHq — dHi for the region Hq. 

We do not have to specify the functional forms of hi and /12 for our theory, but 
assume that the magnitude of hi an d hj is very small. Technically speaking, hi 



and /i2 are nearly flat in the sense of IShimodairal ( 20081 ). Introducing an artificial 
parameter A, a function h is called nearly flat when sup^g^m |/i(e)| = 0(A) and L^- 
norms of h and its Fourier transform are bounded. We develop asymptotic theory 
as A — *■ 0, which is analogous to n — > 00 with the relation A = l/\/n. 

The whole parameter space is partitioned into two regions as HqU Hi = M'"+^ 
or three regions as Hq U HiU H2 = M."^^^. These partitions are treated as disjoint 
in this paper by ignoring measure-zero sets such as Hq n Hi = dHi. Bootstrap 
methods for computing frequentist confidence measures are well developed in the 
literature as reviewed in Section [3l The main contribution of our paper is then 
given in Section |4] for the case of three regions. In Section [5l this new comp utation 



method is used also for Bayesian measures of lEfron and Tibshiranil ( 19981 ). Note 



that the flat prior 7r(/i) = const in the previous section was in fact carefully chosen 
so that 7r(iJ^|y) = p{Hl,\y) for ^. This same 7r(^) led to 7r(iJo|2;) ^ p{Ho\y) for ^. 
Our definition of Hq given in ([9]) is a simplest formulation, yet with a reasonable 
generality for applications, to observe such an interesting difference between the 
two confidence measures. 

Multiscale bootstrap computation of the confidence measures for the three re- 
gions case is described in Section |6l Simulation study and some discussions are 
given in Section [7] and [HI respectively. Mathematical proofs are mostly given in 
Appendix. 

3. Frequentist measures of confidence for testing two regions 
In this section, we review the multiscale bootstrap of IShimodairal (2003) for 



computing a frequentist p- value of "one-sided" test of Hq. Let z = — $~^(a) be the 
inverse function of a = <&(— z). The bootstrap z-value of Hq, defined as z^-^ (-f^o I?/) = 
—^~^{a„2{HQ\y)), is convenient to work with. By multiply ing a to it, az„2(H n\y) 
is called the normalized bootstrap z- value. Theorem 1 of IShimodairal (J2008h . as 



reproduced below, states that the z-value oi p{HQ\y) is obtained by extrapolating 
the normalized bootstrap z-value to ct^ = — 1, or equivalently n' = —n in terms of 
X*. 

Theorem 1. Let Hq be a region of mU\) with nearly flat hi. Given Hq and y, 
consider the normalized bootstrap z -value as a function of a^ ; We denote it by 
■0((7^) — ~a^~^{acr2{HQ\y)). Let us define a frequentist p-value as 

(11) p(H^|2/) = $(_v;(-l)). 
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Figure 1. Region Hq C M™"*"^ is the shaded area between surfaces 
dHi and dH2. 



and assume that the right hand side exists. Then for fj, G OHq and < a < 1, 

(12) P(pK|y)<a|/i)=a + 0(A3), 

meaning that the coverage error, i.e., the difference between the rejection probability 
and a, vanishes asymptotically as X —f 0, and that the p-value, or the associated 
hypothesis testing, is "similar on the boundary" asymptotically. 

Proof. Here we show only an outhnc of the proof by allowing the coverage error of 
O (X^), instead of_0(A^), in p^ . This is a brief summary of the argument given 
in Shi modairal ( 20081 ). First define the expectation operator £„2 for a nearly flat 
function h as 

i£,2h)iu) := E,2ihiU*)\u), 
where E„2[-) on the right hand side denotes the expectation with respect to (O, 
that is, for Y* = {U* ,V*) with 

U*\u^Nr^{u,a''l^), V*\v^N{v,a^). 
Using the expectation operator, we next define two quantities 

V + E„2hi{u) hi{U*) — £^2hi{u) 

a G 

and work on the bootstrap probability as 

a„2{%\y) = P,2(V* < -h^{U*)\v) 
= E„2 mzi + ei)\u) 

= £;,2($(zi) + 0(zi)ei|w) + O(A') 

(13) = $(zi) + 0(A2). 

The third equation is obtained by the Taylor series around zi , and the last equation 
is obtained by E^2{ei\u) = 0. Rearranging (fT5|) . we then get the scaling-law of the 
normalized bootstrap z- value as 

(14) iP{a^)=v + £,2hi{u)+0{\^). 

On the other hand, eq. (5.10) of IShimodairal ( 20081 ) shows, by utilizing Fourier 
transforms of surfaces, that (fT^ holds with coverage error O(A^) for a p-value 
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defined as 

(15) p{H',\y) = M-v - £-ih,{u)) + 0(A2). 

The proof completes by combining p^ and P^ . D 

A hypothesis testing is to reject Hq when observing p{HQ\y) < a for a specified 
significance level, say, a = 0.05, and otherwise not to reject Hq. The left hand side 
of (fT2|) is the rejection probability, which should he < a ioi ji £ H'q and > a for 
fi ^ Hq to claim the unbiasedness of the test. On the other hand, the test is called 
similar on the boundary when the rejection probability is equal to a for /x S OHq. 
In this paper, we implicitly assume that p{HQ\y) is decreasing as y moves away 
from Hq. The rejection probability increases continuously as fi moves away from 
Hq. This assumption is justified when A is sufficiently small so that the behavior 
of p(ifoly) is not very different from that for ([2]). Therefore, (|12p implies that the 
p- value is approximately unbiased asymptotically as A ^ 0. 

We can think of a procedure for calculating p{H'f^\y) based on (fTTjl . In the 
procedure, the functional form of i/)(o-^) should be estimated from the observed a„2 's 
using parametric models. Then an approximately unbiased p- value is computed by 
extrapolating '4>{u'^) to ct^ = — 1. Our procedure works fine for the particular Hq 
of ©, because ip{a'^) = y^+i and piH^^ly) = $(-y„+i) = $(-V'(-l)). Our 
procedure works fine also for any Hq of (fTU]) when the boundary surface OHq is 
smooth. The model is given as V'(^^) = Po + Pia"^ + /32cr^ + P^a^ + ■ ■ ■ using 
parameters cp = (/3o,/?i, ■■ ■), and thus an approximately unbiased p- value can be 
computed by p{HQ\y) = $(— /3o +/3i — /32 +/53 — ■ • • )■ It may be interesting to know 
that the parameters are interpreted as geometric quantities; /3o is the distance from 
y to the surface BH'q, j3i is the mean curvature of the surface, and (3j, j > 2, is 
related to 2j-th derivatives of hi. 

However, the series expansion above does not converge, i.e., '0(— 1) does not 
exist, when BHq is nonsmooth. For example, ipia"^) = (3q -\- (3i\a^ serves as a good 
approximating model for cone-shaped H'q^ for which '0(— 1) does not take a value 
of M. This observation agrees with the fact that an unb iased test does not exist for 



cone-shaped Hq as indicat ed in the ar g umen t of lLehmanni (|1952i ). Instead of (fTTj) . 



the modified procedure of lShimodairal (|2008f ) calculates a p- value defined as 



for an integer fc > and a real number ctq > 0. This is to extrapolate V'(o'^) 
back to CT^ = — 1 by using the first k terms of the Taylor series around Gq. The 
coverage error in p^ should reduce as k increases, but then the rejection region 



viol ates the desired property calle d monotonicity in the sense of iLehmannI (|1952[ ) 



rty 

m 



and iPerlman and Wul (|l999l . l2003f ). For taking the balance, we chose fc = 3 and 



(Tq = 1 for numerical examples in this paper. 

4. Frequentist measures of confidence for testing three regions 

The following theorem is our main result for computing a frequentist p-value of 
"two-sided" test of Hq. The proof is given in Appendix lA.il 

Theorem 2. Let Hq he a region of (Qi with nearly flat hi and /i2. Given Hq and y, 
consider the approximately unbiased p-value p{Hi\y) by applying Theorem[l\to Hi 
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for i = 1,2. Assuming these two p-values exist, let us define a frequentist p-value 
of Hq as 

(17) p{Ho\y) = 1 ~ \p{H^\y) ~ p{H2\y)\. 

For example, |J7p holds for the exact p-value of Q) defined in Section[li Then for 
11 e dHo = dHi U dH2 andQ<a<l, 

(18) P{p{H^\Y)<a\^x)=a + 0{\^), 

meaning that p{Hf)\y) is approximately unbiased asymptotically as A ^ 0. 

For illustrating: the meth odology, let us work on the spherical shell example of 
lEfron and Tibshiranil ( 1998[ ). for which we can still compute the exact p-values to 



verify our methods. The region of interest is H^ : 02 < \\fJ.\\ < ai as shown in 
Panel (a) of Fig. [2l We consider the case, say, 

m + 1 = 4, ai = 6, 02 = 5, \\y\\ — 5.9, 

so that this region is analogous to ([5]) except for the curvature. The exact p-value 
for Hi : ||/i|| > ai is easily calculated knowing that ||1"||^ is distributed as the chi- 
square distribution with degrees of freedom m + l and nonccntrality ||/i|p. Writing 
this random variable as Xm-i-i(l!MlP)i the exact p-value is p{Hi\y) = PiXm+iio-l) — 
II2/IP) = 0.362, that is, the probability of observing ||y|| < \\y\\ for ||/i|| — ai. 
Similarly, the exact p-value for H2 : \\n\\ < 02 is p{H2\y) = P {xin+ii^-V} ^ II?/IP) = 
0.267. In a similar way as for ([3|), the exact p-value for Hq is computed numerically 
asp(iJoly) = 0.907, although the procedure is a bit complicated as explained below. 
We first consider the critical constants ci and C2 for the rejection regions i?i = {y | 
llyll < oi — ci} and R2 — {y \ \\y\\ > 02 +C2}. By equating the rejection probability 
toafor^e9i/o,thatis, P(x^+i(a|) < {ai-ci)^) + P{xl,+i{af) > (oa + ca)^) = a 
for z = 1, 2, we may get the solution numerically as ci — 1.331 and C2 = 1.903 for 
a — 0.05, say. The p-value is defined as the infimum of a such that Hq can be 
rejected. 

To check if Theorem [5] is ever usable, we first compute p7)) using the exact 
values of p(i7i|y) andp(i/2|y)- Then we get p(i/o|2;) = 1 - (0.362- 0.267) = 0.905, 
which agrees extremely well to the exact p{Ho\y) = 0.907. The spherical shell 
is approximated by ([9|) only locally in a neighborhood of y but not as a whole. 
Nevertheless, Theorem [2] worked fine. 

We next think of the situation that bootstrap probabilities of Hi and H2 are 
available but not their exact p-values. We apply the procedure of Section [3] sepa- 
rately to the two regions for calculating the approximately unbiased p-values. To 
work on the procedure, here we consider a simple model '(/'(cr^) = (3q + /3i cr^ with 
parameters ip — (/3o,/3i) for 

(19) a„2{H[,\y)^^-i;{a')/a). 

Let 'ipi{a'^) be the normalized bootstrap z-value of Hi for i = 1,2. By assuming 
the simple model for %lji{a'^), we fit a„2{Hi\y) = $(— V'i(o'^)/o') to the observed 
multiscale bootstrap probabilities of Hi for estimating the parameters. The actual 
estimation was done using the method described in Section [6.31 but we would like 
to forget the details for the moment. We get /3o = 0.101, /3i = —0.258 for Hi, 
and similarly J3o — 0.889, jSi = 0.286 for 7^2- /3o's are interpreted as the distances 
from y to the boundary surfaces, and the estimates agree well to the exact values 
/3o = 0.1 for Hi and [3q ~ 0.9 for iJa- Then the approximately unbiased p-values 
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Figure 2. (a) Spherical shell region, (b) Cone-shaped region (Section [T]). 



are computed by ^ as p(-ffi|y) = $(-0,101 - 0.258) = 0.360 and p(i?2|2/) = 
$(-0,889+0.286) = 0.273, and thus ^ giyesp{Ho\y) = l-(0.360-0.273) = 0.913, 
which again agrees well to the exact p{HQ\y) = 0.907. 

We finally think of a more practical situation, where the bootstrap probabilities 
are not available for Hi and H2, but only for Hq. This situation is plausible in 
applications where many regions are involved and we ar e not sure which o f them 
can be treated as Hi or H2 in a neighborhood of y; See lEfron et all ( 19961 ) for an 



illustration. We consider a simple model 'i/'i(c^) = (3o+/3ia^, ^2(0"^) = d—fio — fiia^ 
with parameters (^ = (/3o, /?i, d) for 

(20) a„^{Ho\y) = I - [H-M<y^) I <y) + H-M<y^) / <y)) 

by assuming that the two surfaces are curved in the same direction with the same 
magnitude of curvature \(3i\. For estimating ip, ([20|) is fitted to the observed mul- 
tiscale bootstrap probabilities of Hq with constraints /3o > —d/2 and d > 0, and 
(p is obtained as /3o = 0.089, /3i = —0.199, d = 0.995. Then the approximately 
unbiased p-values are computed by dni as p{Hi\y) = $(-0,089 - 0.199) = 0.387 
and piH2\y) = $(-0,995 + 0.089 + 0.199) = 0.240 and thus ^ gives p{Ho\y) = 
1 - (0.387 - 0.240) = 0.853. This is not very close to the exact p(i?o|y) = 0.907, 
partly because the model is too simple. However, it is a great improvement over 
ai{Ho\y) = P{a\ < x^+i(ll2/f ) < «i) = 0.320. 

5. Bayesian measures of confidence 

Choosing a good prior density is essential for Bayesian inference. We consider 
a version of noninformative prior for making the posterior probability acquire fre- 
quentist properties. 

First note that the sum of bootstrap probabilities of disjoint partitions of the 
whole parameter space is always 1. For the two regions case, 0:0.2 {Hl^\y)+a„2 [Hi \y) = 
1, and thus az^^iH^^ly) + az„2{Hi\y) = 0. Therefore p(i/o|y) +p(ffi|2/) = 1 for the 
approximately unbiased p- values computed by (fTTj) . suggesti ng that we may think 
of a p rior so that p[HQ\y) = Tr{HQ\y). This was the idea of lEfron and Tibshirani 



(J1998I ) to define a Bayesian measure of confidence of Hq- Since each of Hi and H 



2 



can be treated as Hq by changing the coordinates, we may assume a prior satisfying 
(21) TriH,\y)^piH,\y), i = l,2. 
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It follows from J2i=o '^(^ilv) = 1 t^^^^ 

(22) 7r{Ho\ij) - 1 - {p{H,\y)+p{H2\y)). 

Priors satisfying (|2ip are call ed probability matching priors. The theory has been 



rriors satisiymg l|zl|) are call ed probability matcnmg priors. 1 ne tneory nas been 
developed in literature (Peerg, Il965t iTibshiranil . Il989t IPatta and Mukerieei 120041 ) 



for posterior quantiles of a single parameter of interest. The examples are the flat 
prior 7r(/z) = const for the flat boundary case in Section [1] and 7r(^) ex H/iH"™ for 
the spherical shell case in Section ID 

Our multiscale bootstrap method provides a new computation to 7r(_ffo|2/). We 
may simply compute ([22]) with the p{Hi \y) and p{H2\y) used for computing p{HQ\y) 
of (|17p . Although we implicitly assumed the matching prior, we do not have to 
know the functional form of 7r(/i). For the spherical shell example, we may use the 
exact p{Hi\y) and p{H2\y) to get p{Ho\y) = 1 - (0.362 + 0.267) = 0.371, or more 
practically, use only bootstrap probabihties of Hq to get p{Ho\y) = 1 — (0.387 + 
0.240) = 0.373. 

6. Estimating parametric models for the scaling-law of bootstrap 

probabilities 

6. 1. One-step multiscale bootstrap. We first recall the estimation procedure 



of lShimodairal (J2002L l2008l ) before describing our new proposals for improving the 
estimation accuracy in the following sections. 

Let f{a'^\(p) be a parametric model of bootstrap probability such as P^ for 
Hq or (|20p for Hq. As already mentioned in Section [l] the model is fitted to the 
observed Ci/Bi, i = 1, . . . , M. Since Ci is distributed as binomial with probability 
/(crf|(^) and Bi trials, the log-likehhood function is ^((/?) — X^i^il^'i log/(o'f |(/j) + 
{Bi — Ci)log(l — /(cr||<p))}. The maximum likelihood estimate ip is computed 
numerically for each model. Let dmnp denote the number of parameters. Then 
AIC — ~2t{(p) + 2 dim (/J may be compared for selecting a best model among 
several candidate models. 



6.2. Two-step multiscale bootstrap. IShimodairal ( 2004 ) has devised the multistep- 



multiscale bootstrap as a generalization of the multiscale bootstrap. The usual 
multiscale bootstrap is a special case called as the one-step multiscale bootstrap. 
Our new proposal here is to utilize the two-step multiscale bootstrap for improving 
the estimation accuracy of (/3, although the two-step method was originally used for 
replacing the normal model of ([1]) with the exponential family of distributions. 

Recalling that X* is obtained by resampling from X , we may resample again 
from X* , instead of A", to get a bootstrap sample of size n", and denote it as 
X** = {x**, . . . ,a;** }. We formally consider the parametric bootstrap 

Y**\y* ^ N^+i{y\{T^ ~ <j^)I^+^), 



where r is a new scale defined by t"^ — a"^ = n/n". In IShimodairal ( 20041 ) . only 
the marginal distribution Y**\y ~ A'^m_|_i(j/,T^/,„+i) is considered to detect the 
nonnormality. For the second step, P„^^t2{Y** g Ha\y) = ar2{Ha\y) should have 
the same functional form as P^'i.j.ilY* E Holy) = au2{Ho\y) for the normal model. 
Here we also consider the joint distribution of (^Y*,Y**) given y. It is 2m + 2- 
dimensional multivariate normal with Cov{Y* ,Y**\y) = a^Im+i- We denote the 
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probability and the expectation by Pa-2.^2{-\y) and £'0-2 ,-2 (-ly), respectively. Then, 
the joint bootstrap probability is defined as 

a,2,,2(iJo|y) = P.2,,2(r* eHoA Y** e Ho\y). 

Let (7(cr^, r^|(^) be a parametric model of q;o.2 7.2 (i?Q|j/) or a^2.,.2{HQ\y). To work 
on specific forms of g{<j'^,T'^\ip), we need some notations. Let {X',X") be dis- 
tributed as bivariate normal with mean (0,0), variance V{X') = V{X") = 1, and 
covariance Cov{X' ,X") = p. The distribution function is denoted as $p(ai,6i) = 
P{X' < fli A X" < bi), where the joint density is explicitly given as (pp{ai,bi) — 
(1 - p2)-i/2^((i _ p2)-i/2(^^ _ pai))^{ai). We also define $p(ai, 61; 02, 62) = 
P{a2 <X' <aihh2< X" < bi) = $p(ai, 61) - $p(a2, 61) - $p(ai, 62) + $(02, 62). 
Then a generalization of (fTil) is given as follows. The proof is in Appendix lA. 21 

Lemma 1. For sufficiently small \, the joint bootstrap probabilities for H'q and Hq 
are expressed asymptotically as 

(23) a^2^,2{H^\y) = *p(zi, w;i) + ©(A^), 

(24) a„2^^2{Ho\y) = <Pp{zi,wuZ2,W2) + 0{X^), 

where z\ — '-(v^Ea2h\{u)) j a , w\ = —{v+Er2hi{u))/T, Z2 = ~{v+d+Ecr2h2{u)) / a , 
W2 — —{v + d + £^2h2{u))/T, and p — a/r. 



Thus g{a'^,T'^\ip) is specified for Hq as ([23|) with zi = —ijj{a'^)/a, wi = ~'iP{t'^)/t 
using the ip function of (I19|) . Similarly, g{a'^,T'^\ip) is specified for Hq as (|24|) with 
zi = ■0i(cr2)/cr, wi = iPi{t'^)/t, Z2 = -V'2(o-^)/cr, ^3 = -iP2{t'^)/t using V'l and ip2 
functions of pl)l . 

We may specify M sets of (ct, r), denoted as (ui, ti ),..., (cta/, tm). In our 
numerical examples, (Ti,...,(Ti3 are specified as mentioned in Section [1] and r^'s 
are specified so that rf —erf = 1 holds always, meaning n" = n. For each 
i = 1,...,M, we generate (y*^y**) with (cri,ri) many times, say Bi = 10000, 
and observe the frequencies Ci == #{Y* E Hq), Di = #(y** e Hq), and Ei = 
#(F* e Hq A Y** e iJo)- Note that only one Y** is generated from each Y* here, 
whereas thousands of F**'s may be generated from each Y* in the double bootstrap 
method. The log-hkelihood function becomes i{ip) — '^i^i{Ei\ogg{a'j ,Tf\{p) + 
ia - E.,) \og{f{af\^) - g{al rfW)) + (A - E,) \og{f{T^W) - giaf,ri\ip)) + (B. - 
d - Di + Ei) log(l - f{cTf\(p) - /(rf |(/3) + g{af, r^|(p))}. In fact, we have used this 
two-step multiscale bootstrap, instead of the one-step method, in all the numerical 
examples. 

The one-step method had difhculty in distinguishing Hq with very small d from 
Hq with moderate d but heavily curved dHi. The two-step method avoids this 
idcntifiability issue because a small value of Ei indicates that d is small; It is 
automatically done, of course, by the numerical optimization of £(ip). 

6.3. Higher-order terms of bootstrap probabilities for testing two regions. 

The asymptotic errors of the scaling law of the bootstrap probabilities in (fT3|) and 
([23|) are of order O(A^). As shown in the following lemma, the errors can be reduced 
to O(A^) by introducing correction terms of O(A^) for improving the parametric 
model g{a'^,T'^\(p) of Hq. The proof is given in Appendix I A. 3 1 
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Lemma 2. For sufficiently small \, the bootstrap probabilities for H'q are expressed 
asymptotically as 

(25) a^-{H'^\y) = $(zi + Azi) + 0(A3) 

(26) aMHo\y) = Hwi + /^wi) + 0{X') 

(27) a„2^^2{H^\y) = <^p+Ap{zi + Azi,wi + Awi) + 0{\^), 

where Zi, Wi, and p are those defined in Lemma\^ and the higher order correction 
terms are defined as Azi = ~^ziE„2 .^^[eWu) , Awi = —^wiE„'2.t.2[5Wu), and 
Ap = — i (^pE„2.^2[e\\u) + pE^2 .^2 {Si\u) — 2E^2.^2{eiSi\u)') using 

h.jU*) - £,2h,{u) h,{U**) - £,2hiiu) 

(ZH) ei = , oi — . 

a T 

For deriving a very simple model for Ap, we think of a situation h{u) = {A/y/m)\\u\\ 
{B/m)\\u\\'^ and = 0, and consider asymptotics as tti — > oo. This formulation is 
only for convenience of derivation. The two values A and B will be specified later 
by looking at the functional form of f{a'^\(p). A straightforward, yet tedious, cal- 
culation (the details are not shown) gives ipi'^'^) = const + Aa + Ba^ + 0{m~^) 
and 

Ap = — ^ (A^p[l -p) + 2B''p{t^ - a^) + 2ABa{l - p^)) + 0{m~^^^). 
2m ' 

This correction term was in fact already used for the simple model '0(c^) == 
/3o + /3i cr^ of the spherical shell example in Section 01 where the parameter was 
actually ^p — {f3o,f3i,m) instead oi ip = (/3o,/3i). We did not change the ip{(J^) for 
adjusting Azi and Awi, meaning that zi + Azi, instead of zi, was modelled as 
—ip{a'^)/a. Comparing the coefficients of ■(/'(o'^)i we get A = and B — /3i, and 
thus Ap = -(^i)2(cr/T)(r2 - a'^)/m. When ^ was fitted to Hi, the estimated 
parameter m = 2.83 was close to the true value ttt. = 3. 

For the numerical example mentioned above, we have also fitted the same model 
but Ap = being fixed. The estimated parameters are /3o = 0.101, /3i = —0.256, 
and the p-value is p{Hi\y) = <i>(— 0.101 — 0.256) = 0.361. These values are not 
much different from those shown in Section [4l However, the AIC value improved 
greatly by the introduction of Ap, and the AIC difference was 96.67, mostly because 
improved fitting for the joint bootstrap probability of (P7)) . My experience suggests 
that consideration of the Ap term is useful for choosing a reasonable model of 

7. Simulation study 

Let us consider a cone-shaped region H^ in M? with the angle at the vertex 
being 27r/10 as shown in Panel (b) of Fig. [2l This cone can be regarded, locally in 
a neighborhood of y with appropriate coordinates, as H^ of ^ when y is close to 
one of the edges but far from the vertex, or as H'^ of (jTU]) when y is close to the 
vertex. In this section, the cone is labelled either by Hq or H'q depending on which 
view we are taking. 

Cones in R^ appear in the problem of multiple comparisons of three elements 
Xq,Xi,X2, say, and Hi corresponds to the hypothesis that the mean of Xj is the 
large st among the three (JDuPreez et aj Il985t IPerlman and Wul . l2003l : IShimodairal . 



20081 ). The angle at the vertex is related to the covariance structure of the elements. 
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Although an unbiased test does not exist for this region, we would like to see how 
our methods work for reducing the coverage error. 

Contour lines of confidence measures, denoted p{y) in general, at the levels 0.05 
and 0.95 are drawn in Fig.[3l The rejection regions of the cone and the complement 
of the cone are R — {y\p{y) < 0.05} and R' = {y\p{y) > 0.95}, respectively, 
at a = 0.05. We observe that p{y) decreases as y moves away from the cone in 
Panels (a), (b), and (c); See Appendix IBJ for the details of computation. On the 
other hand, Figs. S] and [5] show the rejection probability. For an unbiased test, it 
should be 5% for all the /i G OHq so that the coverage error is zero. 

In Panel (a) of Fig.[3l p{y) = ai{Ho\y) is computed by the bootstrap samples of 
cr^ = 1 . This bootstrap probability, labelled as BP in Fig. [H is heavily biased near 
the vertex, and this tendency is enhanced when the angle becomes 27r/20 in Fig. [5l 

In Panel (b) of Fig. [3l p{y) — p{HQ\y) is computed by regarding the cone as Hq 
of (jlOp. The dent of R and the bump of R' become larger than those of Pane l (a) 
of Fig. [3] near the vertex, confirming what we observed in IShimodairal ( 20081 ). As 



seen in Figs. H] and [H the coverage error oip{HQ\y), labelled as "one sided" there, 
is smaller than that of BP. 

In Panel (c) of Fig. [3l p(iJo|y) is also computed by regarding the cone as Hq 
of ^, and then one of p(ifo|j/) and p{Ho\y) is selected as p{y) by comparing the 
AIC values at each y. This p{y), labelled as "two sided Freq" in Figs. S] and [51 
improves greatly on the one-sided p- value. The coverage error is almost zero except 
for small ||/i||'s, verifying what we attempted in this paper. The corresponding 
Bayesian posterior probability, labelled as "two sided Bayes," performs similarly. 
Note that the coverage error was further reduced near the vertex by setting simply 
ply) = plHoly) without the model selection (the result is not shown here); However, 
the shapes of R and R' became rather weird then in the sense mentioned at the 
last paragraph of Section [3l 

8. Concluding Remarks 

In this paper, we have discussed frequentist and Bayesian measures of confidence 
for the three regions case, and have proposed a new computation method using the 
multiscale bootstrap technique. In this method, AIC played an important role for 
choosing appropriate parametric models of the scaling-law of bootstrap probability. 
Simulation study showed that the proposed frequentist measure performs better for 
controlling the coverage error than the previously proposed multiscale bootstrap 
designed only for the two regions case. 

A generalization of the confidence measures gives a frequentist interpretation 
of the Bayesian posterior probability as follows. Let us consider the situation of 
Theorem [51 If we strongly believe that fi ^ H2, we could use the one-sided p- value 
p{HQ\y) — 1 — p{Hi\y), instead of the two sided p(iJo|2/)- Similarly, we might use 
l—p{H2\y) if we believe that n ^ Hi. By making the choice "adaptively," someone 
may want to use p^^' (Holy) = l^raax{p{Hi\y),p{H2\y)), although it is not justified 
in terms of coverage error. By connecting p^^' {Ho\y) and p{Ho\y) linearly using an 
index s for the number of "sides," we get 

p'^'HHoly) ^ niHoly) + srmn{piHi\y),p{H2\y)). 

It is easily verified that p{Ho\y) = p^^\Ho\y) and 7r{Ho\y) = p(°)(_ffo|y), indicating 
that the Bayesian posterior probability defined in Section [5l can be interpreted. 
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Figure 3. Contour lines p{y) — 0.05 and p{y) — 0.95. The cone- 
shaped region Hq is rotated so that one of the edges is placed along 
the X-axis. Solid curves are drawn for (a) the bootstrap probability 
with 0-^ = 1, and for (b) the frequentist p- value for "one-sided" 
test. In Panel (c), p{y) is switched to the frequentist p- value for 
"two-sided" test when appropriate. The dotted curve in Panel (c) 
is for the Bayesian posterior probability. 
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Figure 4. (a) Rejection probability of the cone, and (b) that of 
the complement of the cone. The angle at the vertex is 27r/10. 
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Figure 5. (a) Rejection probability of the cone, and (b) that of 
the complement of the cone. The angle at the vertex is 27r/20. 



interestingly, as a frequentist p- value of "zero-sided" test of Hq- Although we 
have no further consideration, this kind of argument might lead to yet another 
compromise between frequentist and Bayesian. 

Our formulation is rather restrictive. We have considered only the three regions 
case by introducing the surface /12 in addition to the surface hi of the two regions 
case. Also these two surfaces are assumed to be nearly parallel to each other. It 
is worth to elaborate on generalizations of this formulation in future work, but 
too much of complication may result in unstable computation for estimating the 
scaling-law of bootstrap probability. AIC will be useful again in such a situation. 

Appendix A. Proofs 

A.l. Proof of Theorem [2j First we consider rejection regions of testing Hq for 
a specified a by modifying the two rejection regions of ([3]). Since hi and /12 are 
nearly flat, the modified regions should be expressed as Ri — {{u,v) \ v > c — 
ri{u), u G IR™} and R2 — {{u,v) \ v < —d — c — r2{u), u € M™} using nearly 
flat functions ri and r2- The constant c is the same one as defined in @. Write 
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a = 0(c), b = (t>{c + d) for brevity sake. We evaluate the rejection probability for 
IX e dHi U dH2- Let fi e dHi for a moment, and put /i = {6, —hi{6)). By applying 
the argument of P^ to i?i but ([7]) is replaced by (IT|), we get P{Y E -Ri|m) = 
1 - $(c- £iri (9) + hi (9)) + 0(A2) = $(-c) + a(£iri (6*) - hi (9)) + 0{\^). The same 
argument applied to i?2 gives P{Y e i?2|Ai) = ^{-d~c-£ir2{9) + hi{9j) + 0{X'^) == 
$(— d— c) + b{—£ir2{9) + hi{9)) + O(A^). Rearranging these two formula with the 
identity 

(29) P{Y £ Rilfi) + P{Y e R2\fi) -- a 

for an unbiased test, we get an equation a{£iri{9) — hi{9)) + b{~£ir2{9) + hi{9)) = 
O(A^). By exchanging the roles of n and r2, the equation becomes b{£iri{9) — 
h2{9)) + a{-£ir2{9) + h2{9)) = 0{\^) ioi fi e dH2 with fj. = (6*, -d- /i2(6')). These 
two equations are expressed as 

For solving this equation with respect to ri and r2, first apply the inverse matrix 
of the 2x2 matrix from the left in ((30)) . and then apply the inverse operator of £i 
so that 

(31) ( ^-f] ]=^(l M f f-''!f\ )+0{X'). 
^ ' \ T2{u) ) a + b\ b a J \ £-ih2{u) J ^ ' 

Next we obtain an expression of p- value corresponding to the rejection regions. 
p(ifo|y) is defined as the value of a for which either of y S dRi and y £ 9i?2 holds. 
Note that ri, r2, and c depend on a. Let us assume y G dRi and thus c = v + ri(u) 
for a moment. Write a' = (j){v) = a + 0{X), b' = (t){v+d) = b+0{X) for brevity sake. 
Recalling H]), p{Ho\y) = $(-c) + $(-d-c) = <^{-v - ri{u)) + <^{-d- v - ri(u)) = 
$(— u) + <i>(— d— v) — [a' + 6')ri(w) + O(A^), where ri(w) in ([3T|l can be expressed 
as 

'^i(") = ^XTT^-i'^i^^) + ^^^-i'*2(«) + 0{X^). 
a' + b' a + b' 

Therefore, p{Ho\y) = $(-«) + $(-d - u) - a'£_i/ii(M) - b'£^ih2(u) + 0{X^) = 

$(-w - £^ihi{u)) + $(-d - V - £-ih2{u)) + 0{X^). By applying ^ to i?i 

and H2^ respectively, we get p{Hi\y) — $(« + £^ihi{u)) + 0{X?) and p(iJ2|2/) = 

$(-t.-d-£_i/i2(w))+0(A2), and thus p(7Io|2/) = l-p{Hi\y)+p{H2\y)+0{X^). By 

exchanging the roles of -ffi andiJ2, we havep(-ffo|j/) = l—p{H2\y)+p{Hi\y)+0{X^) 

for y £ 9i?2- By taking the minimum of these two expressions of p(iJo|y), we finally 

obtain ([TT]). This p-value satisfies ^^ with error O(A^), and thus ([11]) holds. 

A. 2. Proof of Lemma [TJ The argument is very similar to (fT5|) in the proof of 
Theorem[T] Given v, u*,u**, the joint distribution of X' = {V* ~ v)/a and X" — 
{y** -v)It is $p. Therefore, Pa^r-^^V* < -hi{u*)AV** < -hi{u**)\v,u*,u**) = 
P„2^^2{X' < zi + ei A X" < wi + 5i\v,u*,u**) = ^p{zi + ei, wi + (5i), where ei 
and 5i are defined in (p8| . Taking the expectation with respect to (t/*, [/**), we 
have a,2^^2{H'o\y) = P,2^^2{V* < -hi{U*) A V** < -hi{U**)\y) = E,2^^2{^p{zi + 
eijWi -\-5i)\u). For proving ([23| . considering the Taylor series around (2:i,wi), we 
obtain 

(32) E,2^,2 Up{zi,wi) + ^£1 + 1^-51 h] + 0(A2) 
with E„2.^2{ei\u) ~ E„2^^2[5i\u) = for completing the proof. 
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Next we show ([M)) . The conditional probabihty given v,u*,u** is Pg.2.,-2{—d — 
h2{u*) <V* < ~hi{u*)A-d-h2iu**) < V** < ~hi{u**)\v,u*,u**) = P^2^^2{z2 + 
62 < X' < zi + ei A W2 + $2 < X" < wi +6i\v,u*,u**) = $p(zi + ei,wi + Si; Z2 + 
e2,W2 + 82), where 

h2iU*) - E,2h2(u) ^ h2iU**) - Er2h2{u) 

£2 = , 02 = . 

(7 T 

Taking the expectation with respect to {U* , U**), weha,veacr2.^2{Ho\y) = Pa^^r^i—d- 
h2{U*) <V* < -hi{U*) A -d-h2{U**) < V** < -hiiU**)\y) = E^2^^2{^p{zi + 
ei,wi + Si; Z2 + £2, W2 + S2)\u). We only have to consider the Taylor series 

/ 9$ 9$ 9$ 9$ \ 

E„2 2 <^ {zi,wi; Z2,W2) + ^ei + ^61 + —-^€2 + T-^62\u] +0{X^) 
\ ozi awi 0Z2 0W2 J 

with E„2.,.2{f.i\u) = E„2,^2{6i\u) = 0, i = 1,2 for completing the proof. 

A. 3. Proof of Lemma [2l By considering a higher-order term of the Taylor series 
in dni), we obtain a^2{H^\y) = E^2{^{zi) + (j){zi)ei - (j){zi)zie\/2\u) + O(A^) = 
$(zi) + 0(zi)Azi + 0{\^) = $(21 + Azi) + 0(A3), proving ^ as well as (^5)1 . 
On the other hand, (j27p is shown by considering higher-order terms of the Taylor 
series in (1321) as 



E„2^^2i^^p{zi,wi) + a/^1+9^'^1 
The proof completes by rearranging the above formula with 

-^-^ = -Zl^ p0p(zi,Wl), =(l)p{zi,Wi), —- =(l)p{zi,Wi). 

ozf ozi oziowi op 

Appendix B. Simulation Details 

The contour lines in Fig. |31 are drawn by computing p- values at all grid points 
(300 X 180) of step size 0.05 in the rectangle area; This huge computation was 
made possible by parallel processing using up to 700 cpus. The computation takes 
a few minutes per each grid point per cpu. Ou r algorithm i s imp lemented as an 
experimental version of the scaleboot package of [Shimodaira| ( 20061 ). which will be 



included soon in the release version available from GRAN. 

The rejection probabilities in Figs. IH and IS are computed by generating y ac- 
cording to (m for 10000 times, and then counting how many times p{y) < 0.05 
or p{y) > 0.95 is observed. This computation is done for each ^ g OHq with the 
distance from the vertex ||^|| = 0, 1, . . . , 16, i.e., ^ = (0, 0), (1, 0), . . . , (16, 0) in the 
coordinates of Fig. I3| 

For computing p(iJg|j/) and p{HQ\y) , the two-step multiscale bootstrap described 
in Section W^ was performed with the M = 13 sets of scales {ai,Ti), i = 1, . . . , 13, 
specified there. The parametric bootstrap, instead of the resampling, was used for 
the simulation. The number of bootstrap samples has increased to Bi = 10^ for 
making the contour lines smoother, while it was Bi = 1 0^ in the other res ults. 

For p{HQ\y), we have considered the singular model of lShimodairal (2008!) defined 



as tp{a'^) = /3o+/3i/(l+/32(o'— 1)) for cones, and performed the model fitting method 
described in Section lOl From the Taylor series of this ^(c^) around c = 1, we get 
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A == /3i/32 (3 — 2/?2) , B = (3i{f32 — l)^ for computing the higher order correction term 
Ap. We have also considered submodels by restricting some oiip — (/3o, /?i, /32, "^) to 
specified values, and the minimum AIC model is chosen at each y. The frequentist 
p- value is computed by (|16p with fc = 3 and CTq = 1. 

For p^HqIh), we have considered the same singular model for the two surfaces 
by assuming they are curved in the opposite directions with the same magnitude 
of curvature. More specifically, the two tp functions in ([20|) are defined as ■(/'i(cr^) — 
/3o+/3i/(1 + /32(ct-1)) andf/'2(cr^) = d~/3o + /3i/(l + /32(cr-l)). The parameters </? = 
{Po, f3i,l32,d) are estimated by the model fitting method described in Section W^ 
Submodels are also considered and model selection is performed using AIC. The 
frequentist p- value is computed by ()17p . and the Bayesian posterior probability is 
computed by ([22|) . 

The rejection probabilities o f other two c ommo nly used measures are shown 
only for reference purposes; See IShimodairal ( 20081 ) for the details. The rejection 



probability of the multiple comparisons, denoted MC here, is always below 5% in 
Panel (a), and the coverage error becomes zero at the vertex. On the other hand, 
the rejection probability of the z-test is always below 5% in Panel (b), and the 
coverage error reduces to zero as ||/i|| — > cx). 
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